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Contact between an elastic manifold and a rigid substrate with a self-afiine fractal surface is 
reinvestigated with Green's function molecular dynamics. The Fourier transforms of the stress and 
contact autocorrelation functions (ACFs) are found to decrease as |q|~'' where q is the wavevector. 
Upper and lower bounds on the ratio of the two correlation functions are used to argue that they 
have the same scaling exponent /i. Analysis of numerical results gives ^ = 1 -\- H, where H is 
the Hurst roughness exponent. This is consistent with Persson's contact mechanics theory, while 
asperity models give = 2(1 + H). The effect of increasing the range of interactions from a hard 
sphere repulsion to exponential decay is analyzed. Results for exponential interactions are accurately 
described by recent systematic corrections to Persson's theory. The relation between the area of 
simply connected contact patches and the normal force is also studied. Below a threshold size the 
contact area and force are consistent with Hertzian contact mechanics, while area and force are 
linearly related in larger contact patches. 
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I. INTRODUCTION 



The distribution of pressures in the contact between two solids has a crucial influence on the amount of friction 
and wear that occur when the solids slide against each other. As a consequence, there is great interest in reliable 
predictions for the dependence of the distribution on the externally imposed load L, the mechanical properties of 
the solids in contact, and their surface topographies^. Of course, it would be desirable to have theories at hand that 
allow one to calculate these distributions analytically or numerically at a moderate amount of computing effort. While 
solving numerically the full elastic or plasto-elastic behavior of contacting solids has become an alternative to studying 
analytical theories, numerically exact approaches incorporating long-range elastic deformation remain challenging to 
pursue. — This is because treating the wide range of length scales present in most surface topographies places great 
demands on computing time and memory. 

Traditionally, many contact mechanics predictions were based on models following Greenwood and Williamson's 
(GW) seminal paper^, in which the contact between two solids is treated as the sum of non-interacting, single-asperity 
contacts. The crucial, geometric properties entering such theories are the statistics for asperity height and curvature. 
One of the most detailed treatments is by Bush et ali^ who included a distribution of curvatures and elliptical 
asperities. In the last decade, Persson and collaborators have pursued a different approach'*'^'^'^, in which the height 
autocorrelation function (ACF) is the geometric property that determines the pressure distribution and thereby the 
area of contact. 

A central quantity in both GW and Persson's theory is the ratio of the real area of microscopic contact to the 
macroscopic projected area of the surfaces Aq. In the limit of low loads, both theories find that 
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where k is a dimensionless constant of proportionality, ctq = L/Aq is the mean pressure normal to the interface, 
\/ (S/h?) is the root mean-square gradient of the height profile h, and the effective elastic modulus E' — E/{1 — v'^), 
where E is the Young's modulus a nd v the Poisson ratio. The generalization of GW by Bush et al^^, predicts n — \/27r, 
while Persson theory yields k = ^/sJtt. Numerically exact calculations for the relative contact area Ac/Aq for solids 
with a self-affine surface topography are, give or take, half way between the two theoretical prediction a^'^'^" . 

Persson's theory also provides a prediction for the functional dependence of the stress distribution P((t), which can 
be described as the sum of a Gaussian of width Act centered at the macroscopic pressure ao and a mirror Gaussian 
centered at —ctq: 
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A 5-function contribution is added to P(cr) so that the integral over P(cr) from cr = to a = oo yields unity. The 
prefactor to this J-function contribution (divided by the normal macroscopic stress) can be interpreted as the relative 
area of the projected surface that is not in contact with the counterface. For the relevant non-negative pressures, the 
superposition leads to P(cr) oc a for small a and Gaussian tails at large a. GW theory gives distributions with similar 
limiting behavior. Numerical solutions for P(cr) are qualitatively consistent with Eq. [2]when the parabolic peaks of 
asperities are resolved^i^^, but some quantitative discrepancies remain as illustrated in the next section. 

GW-type theories and Persson theory give quite distinct predictions for the spatial distribution of contact and 
pressure. Contact in GW is based on overlap of undeformed surfaces. This allows one to relate the contact ACF, 
Cc(Ar), to the surface topography. Cc(Ar) gives the probability of finding contact at a coordinate r if there is contact 
at a coordinate r' = r + Ar that is located at a distance Ar — |Ar| away. In contrast, Persson's theory contains a 
prediction for the stress ACF, Co-(Ar'). Both GW and Persson theory predict power law scaling for the correlation 
functions, but with very different exponents. In a numerically exact calculation, Hyun and Robbins found that Cc 
and Co- appear to decrease algebraically with Ar in such a way that the two functions are essentially proportional 
to one another—. The values for the exponents that describe the decay of the autocorrelation functions were again 
found to be, give or take, half way between the theories. However, due to computational limitations, the uncertainty 
on the exponents was relatively large and only H = 0.5 was studied. 

Here, we would like to reassess these correlation functions with Green's function molecular dynamics (GFMD)ii, 
which allows one to address larger system sizes than with finite-element methods or multiscale approachesiS,. One of 
our main goals in these calculations is to assess whether the relatively precise values of k predicted by GW and Persson 
theory are to a certain degree fortuitous or if the theories also predict other computable observables, specifically stress 
and contact ACFs, with similar accuracy, i.e. of order twenty percent. We derive approximate and exact bounds 
on the ratio of these ACFs and their integrals that supplement the numerical results. The work presented here 
also includes a test of the claim that corrections to Persson's theory can be derived with the help of a systematic 
expansion scheme derived recently by one of us''^'^. The expansion has so far only been worked out (to harmonic 
order) for walls that repel each other with forces that increase exponentially when the distance between the surfaces 
decreases. Therefore, we also include comparison to numerically exact simulations based on exponentially repulsive 
walls. Finally, we conduct an analysis of connected contact patches for hard wall interactions to see if these regions 
show the characteristic behavior of Hertzian contact mechanics, as one would expect according to an overlap theory 
of purely repulsive, corrugated walls. 

The remainder of this paper is organized as follows. In section |TT1 the GFMD is quickly introduced. In section [1111 
bounds on the correlations are derived and numerical results for scaling behavior presented. Conclusions will be drawn 
in section HVl 

II. METHOD 

In this work, we model elastic, frictionless contact between two solids with self-affine surfaces. Use is made of the 
mapping of such a system onto contact between a flat elastically-deformable solid and a rigid, corrugated substratei^. 
The Green's function molecular dynamics (GFMD) methodii was used to solve for the elastic response of the de- 
formable solid. Details of this method are provided in Ref.^^ and we only give the parameters of the model here. 
We chose both Lame constants to be unity, which is equivalent to a Young's modulus oi E — 5/2, a bulk modulus 
oi K = 5/3, and a Poisson ratio oi ly = 1/4. In what follows, most stated quantities will be dimensionless, but we 
take the Lame constants as our unit of pressure. The continuum equations have no intrinsic length scale, so we will 
normalize lengths by the lateral resolution a of the GFMD. 

The surface topography of the rigid substrate was a self-afRne fractal. Surfaces with the desired value of the Hurst 
roughness exponent H were created using a Gaussian filter technique for the Fourier components of the height profile 
h{q)^^. The long wavelength cutoff of fractal scaling Ai is always identical to the length of our periodic simulation cell 
in both lateral directions. The effect of reducing Ai is discussed in Ref?i2.. Unless noted X\/a = 4096. The effective 
depth of the deformable solid is also equal to Ai. The short wavelength cutoff A.; was varied from a to 64a. The 
wavevectors associated with the cutoffs and resolution are denoted as qi = 27r/Ai, Qs = 2tt/Xs, and (/a — 27r/a. The 
magnitude of the Fourier components was adjusted to maintain the same mean-square height gradient (Vft-^) = 0.031 
for all As and H. 

A wide range of parameters was considered. The roughness exponent was varied from 0.3 to 0.8, which covers the 
range of values typically reported in experiments ^^'^^i^^'^^ . The load was varied between 0.001 and 0.256, resulting in 
fractional contact areas /c = Ac/Aq between 0.02 and 0.96. For each case we evaluated the spatial variation of the 
contact pressure (j{r). Then the total contact area Ac, probability distribution P{(j) of local pressures, and contact 
and pressure correlation functions were calculated. Since P(cr) enters the analysis below, we show typical results in 
Figure [TJ along with the analytic expression of Eq. [2l For this case As > > a, implying that the top of each asperity 
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FIG. 1: Pressure distribution P{cr) for a system with Ai — 64As = 4096a at a relative contact area of fc — Ac/Aq = 0.08. The 
line is a fit with Eq. [2] From the area and load we find k = 2.06 (Eq. [l]) with an absolute uncertainty due to finite resolution 
of less than 0.1. 

is resolved into a smooth parabolic peak. As a result, P{(j) drops as a decreases to zero. Extending the roughness to 
X = a changes the form of P{<7) at large and small a, but does not change the power law behavior of the correlation 
functions at small q that are our focus herein. 
The contact ACF, Cc(Ar), is defined as 

C,(Ar) = (e{a(r)}e{a(r')}) (3) 

where Ar = |r — r'| and a{r) is the normal component of the stress at position r.^^ The Heaviside step function, 8(...), 
is zero for negative and unity for positive arguments. However, unlike the usual convention, we choose 6(0) — 0, i.e., 
the step function is unity only at those locations where there is contact. The stress ACF is similarly defined as: 

a(Ar) = (a(r)a(r')) , (4) 

The ACF are most conveniently calculated by Fourier transforming. We choose the normalization so that 

a(g)-(a*(g)a(<z)). (5) 

III. RESULTS 

A. Relation between contact and stress autocorrelation functions 



Traditional models ignore correlations in surface displacement due to elastic deformation so that the location of 
contacting regions is determined solely by the local height. We will refer to such models generically as "overlap 
models" . A particular overlap model is the bearing area model^^ in which contact occurs wherever the undeformed 
solids would overlap. For the case of rough on flat considered in our calculations, this corresponds to the region where 
the height of the rough solid is above a threshold value. In the GW model and extensions^ii^, the same criterion is 
used to determine which asperities are in contact, and the corresponding load is obtained from the force needed to 
remove overlap. Since contact only depends on the local height in such models, it is relatively easy to construct the 
contact morphology for any given surface topography. One can then calculate the relative contact area and Cc(Ar). 

Since the location of contacts is entirely determined by the height in the bearing area and GW model, Hyun and 
Robbins argued that Cc should have the same scaling as the height-height correlation C\i^^: 

Cc{q) ^ C^{q) ^ q-^^^+"y (6) 

Their numerical results for H — 0.5 were consistent with this relation. We are not aware of any specific predictions 
for the stress ACF in the GW model, although it might also follow Ch since the stress on asperities is also determined 
directly by their height. 
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In contrast, Persson's scaling theory does not consider Cc, but contains implicit predictions for the stress ACF. In 
the limit of complete contact, Persson theory becomes exact. The stress correlation function can be determined from 
the height correlation function and the elastic Greens function. The latter scales as g^, yielding Ca ^ q~^^^^ . This 
result is consistent with numerical studies of full contact^-. The original version of Persson theory does not discuss 
stress ACFs in partial contact explicitly. One may thus be tempted to use the full contact approximation for the 
stress ACF at all loads. However, when contact is not complete, only a fraction of the elastic manifold conforms 
to the counterface and so the full contact ACF only provides an upper bound, as the non-conforming parts of the 
manifold do not carry any load. After receiving a preprint of our work, Persson informed us that a correction factor 
needs to be included into the calculation to capture this effeclj^^. This correction factor is the relative contact area at 
a given "magnification", i.e., one needs to correct with the relative contact area /(g) that his theory would predict if 
all roughness for wavevectors of magnitude greater than q were eliminated. (This correction factor had already been 
introduced for the calculation of adhesive interactions in earlier work^.) With this correction factor the exponent 
^„ = \ -\- H in obtained in Persson theory for partial contact. In particular, 

CM^^q^C^il)- fiQ) -q-^'+"^ (7) 




where the first term is the full contact prediction and f{q) ^ q^^^~^^ — 

While the power laws in Equations [6] and [7] are very different, it is not clear how the scaling of the stress ACF and 
contact ACF should be related. In the following we argue that two approximate bounds on the ACF's force them to 
have the same scaling exponents. In particular, 

Cc ( Ar ) < — ( Ar) < ( Ar ) « 2Cc ( Ar) , (8) 

where — {a)^ and (c^)^ are the mean and second moments of the stress averaged over those areas where there is 
contact, i.e., 

jo+ d(jP{a) 

In cases where the stress histogram, P{(j) — {5{a — ^(i")}) can be described by equation ([2]), one can easily find 
that (o'^)c/o'c ^ 2. The same ratio is obtained for the approximately exponential distribution of pressures found for 
surfaces with = . Similar ratios are obtained for the other As considered here, leading us to add the approximate 
equality at the right of Eq. [H 

In the limit Ar — > 0, the stress ACF is exactly equal to the upper limit of Eq. [H Ccr(O) = {(t'^)c- In the large 
Ar limit, the local values of the stress should become decorrelated, and will then equal the lower limit. One 
expects a smooth crossover between the two bounds as Ar increases unless there is a strong correlation at some given 
wavelength. For example, one can construct counterexamples to the bounds such as a perfectly sinusoidal surface 
topography. 

Figure m presents typical numerical results for the ACF's of self-afhne surfaces. Values of Ca lie close to the upper 
bound up to Ar ~ 8a, and then cross over rapidly to the lower bound. A heuristic reason for the more rapid drop in 
Co- than Cc can be constructed as follows. Consider the contribution to these two ACFs that stem from two points r 
and r' that are in the same contact patch, for example, within a single Hertzian contact. The contribution to Cc(Ar) 
will simply be unity, i.e., Cc(Ar) cannot decay within a simply-connected contact region. Conversely, C(j(Ar) can and 
will have a lot of structure, e.g., the correlation between the stress at the center and edge of a Hertzian contact will 
be small. Consequently, a significant fraction of Ccr(Ar) will have decayed on a length-scale Ar that is comparable 
to a typical contact radius, while only a very small fraction of Cc(Ar) will have decayed on that same distance. Note 
that the rapid decay in Fig. [2] starts at Ar/Ag 1/8 which should be comparable to the smallest contacts. Since 
contacts of many different sizes are found, one may conjecture that Ccr(Ar) falls off faster than Cc(Ar) for all Ar. The 
difference should decrease at large Ar as the number of connected clusters with dimension greater than Ar decreases. 

Fig El illustrates that the same upper and lower bounds describe the Fourier transforms of the correlation functions. 
Data for a range of H and relative areas are shown. In each case there is a crossover from the lower bound of Eq. [8] 
at small q (large C) to the upper bound at large q. 

One can derive additional relations for the Fourier transforms of the ACF's that support the bounds quoted above. 
First, the g = values must obey: 

Cc((Z-0) - |Sre(a(r))|2 = [A^AcMo]' (10) 
a(g = 0) = |I]ra(r)|2 = [iVacAcMo]' , (11) 
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FIG. 2: (Color online) Stress autocorrelation function (ACF) Co-(Ar) in real space for H = 0.3 and comparison to the upper and 
lower bound estimates of the stress ACF via equation ^ and the contact ACF Cc. Here Ai = 64As ~ 4096a and Ac/Ao = 0.14. 




C(q) 



FIG. 3: The stress autocorrelation function Ca{q) normalized by as a function of Cc(g). The two straight lines reflect the 
inequalities from equation ((Sj. 



where N is the number of r (or real-space grid points) in the sum. This shows that in the limit q — s- the correlation 
functions satisfy the lower bound in Eq. [8] One can also use a general sum rule over real and reciprocal space to 
write: 



and similarly 



E,a(<z) = S,|a(9)|2 = iVS]r|a(r)|2 = {a^),N^AJAo 



(12) 
(13) 



The last two equations show that the sum over all q of the ACFs is exactly consistent with the upper bound in Eq. 
[51 While this implies that the upper bound must be exceeded for some g, we will see that the sum is dominated by 
large q where the upper bound is nearly obeyed. Our focus is on the power law scaling regime in the opposite limit 
of small q. 
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In order to collapse data for different loads it is useful to recast the above sum rules. It is also helpful to eliminate 
the q = term since we will see that the ACFs diverge in the limit q ^ 0. In addition, the q — term scales as 
{Ac/Aq)"^, while other terms scale linearly with Ac/Aq. Subtracting Eq. [Tl]from Eq. [T3l and rearranging factors we 
find: 

S,^o/cC'e(<z)/(l-/c)Cc(0) = l, (14) 

where /c = Ac/Aq. We will see that scaling Cc in this way removes the dependence on /c at small fc- Evaluating the 
same weighted sum for the stress ACF yields 

S,#o/cC',(<?)/(l - /c)a(O) = [(a')cK' - /c]/[l - fc] > 1. (15) 

Using the same scaling of the two ACF's guarantees that they coincide as g 0. While Ca will decay more slowly 
at large q, the fact that its integral over all q is larger by only a factor of order 2 for any system size implies that Co- 
should be described by the same scaling exponent as Cc. 



B. Comparison to overlap theories 



Figure m compares results for Cc from the full GFMD calculation and the GW model which uses overlap to determine 
contact. Both models yield roughly power law behavior at intermediate q 

Cc(g)ocg-^^ (16) 

However there is a large discrepancy in the numerical values between GFMD (/ic « 1.3) and the bearing area model 
{pc ~ 2.8), which one can conclude from figurelH A similar difference in the values for /Zc is observed for other H and 
As, as reported for H = 0.5 in Refj^^. There is a simple physical reason that overlap models yield larger and thus 
cluster the contact patches too closely. They neglect the fact that an asperity in the vicinity of a very high asperity 
is less likely to come into contact because it is pushed down when that very high asperity comes in contact with the 
counter-surface^iiSi^^. 




FIG. 4: The contact ACF, Cc(g), as a function of wavevector q divided by the long wavelength cutoff qi. Predictions from the 
bearing area model and GFMD calculations are shown for a Hurst roughness exponent of H = 0.3. Here Ai = 2048As — 2048a 
and Ac/Aq = 0.09. The lines reflect power laws in q. 

In all cases, our numerically determined exponent for the bearing area model is consistent with the estimate 
/ic = 2(1 + H), within our uncertainty of ~ 0.25 for this model. In contrast, the GFMD results are consistent with 
fi = 1 + H within an uncertainty of ^ 0.1. The reason for the difference in the uncertainties of our estimates stems 
from the fact that contact geometries associated with large fj,c have large finite-size effects, particularly when fic 
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exceeds twoj21 This is because a large value of /Zc implies significant contributions from long wavelengths, which have 
less sampling than short wavelengths, as can be seen from an analysis of Cc(r): 

Cc(r) = y"d2gC'c(q)expHq-r) (17) 

- / aqq^-f"- / d6lexp(-igrcos(6l)) . (18) 

J a Jo 

For /Xc < 2 the contribution at small q is finite. One can change the integration variable to qr and finds: 

Cc(Ar) - Cc(oo) cx Ar'"'' (19) 

at large Ar with = 2 — /ic. This is the type of behavior we find for our numerical solutions, i.e. Fig. [2l 

The behavior becomes qualitatively different if /ic > 2 as predicted by the /ic = 2(1 + H) > 2 relation for the 
bearing area model. This implies a singular contribution from small q in Eg. 1181 and suggests that Cc(Ar) should 
increase with Ar. The origin of this behavior seems to be related to the distribution of sizes of connected regions in 
the bearing area model. The probability that a cluster will have area Cc is predicted^^ to scale as P(ac) ~ a""^ with 
r = (2 — H/2). Since r < 2, most of the contact area is in the clusters of largest size, and there will be no decay in 
Cc(Ar) on this scale. The dominance of large clusters leads to significant fiuctuations in data for the bearing area 
model that are not seen in the numerical solution with GFMD. We conclude that the bearing area model and thus 
GW provide a very poor description of the contact ACF. 



C. Numerical Determination of /i 



To determine accurate values of the scaling exponents describing the ACF's we maximized the scaling region by 
taking Ag = a. Our results and earlier workl*^ show that resolving the asperities is not important to the large scale 
behavior of interest here. Figure [5] shows results for Ca and Cc sX H — 0.3, H — 0.5 and H = 0.8. In each case, data 
for two different area fractions, corresponding to about 4% and 8% have been collapsed by plotting fcC{q) / (1 — /c)C'(O) 
following Eq. [TH With this rescaling, the data for the two loads are indistinguishable. This is consistent with previous 
numerical work that indicates area fractions of less than 10% exhibit scaling behavior consistent with the asymptotic 
small /c limit^J^iSS. 




FIG. 5: (Color online) Comparison of the stress (open symbols) and contact (closed symbols) ACF's for H = 0.3, 0.5 and 0.8 from 
top to bottom. In each case results for /c near 4% (triangles) and 8% (squares) are collapsed by plotting fcC'{q)/{l — fc)C{0). 
Curves for H = 0.5 and 0.8 are offset vertically downward by successive factors of ten to prevent overlap. Straight lines are fits 
with ^ = 1 + H . In all cases Ai — 2048a and As = a. 
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At small q the stress and contact ACF follow the same scaling. Straight lines show that the data are consistent 
with fic = fJ-a- = ^ + H in this regime. This ansatz for the scaling exponent was motivated by the fact that fi is 
bigger than the value 2H predicted for full contact by an amount that decreases with increasing H (see below). It 
also ensures that /i remains below 2 as H increases to unity, and thus that the real space correlation function has 
nonsingular scaling behavior. If both stress and contact ACF are assumed to have the same exponent in numerical 
fits to the data, then the deviation from 1 + iJ is less than 0.1 (Table However, the fact that Cc always decays 
by about an extra factor of two means that separate fits always give ^.c > jJLa- Given that our scaling range is only 
a decade and a half, the magnitude of this difference is of order log]^Q(2)/1.5 « 0.2. We cannot rule out deviations 
between //o- and /Uc on this scale and it represents the major source of uncertainty in Table HI Given this uncertainty, 
it is possible that may reach or exceed 2 before H reaches unity. As for overlap models, this would imply singular 
contributions from large contacts, and anomalous behavior of the correlation function at large Ar. 





GFMD 


GW 


Persson 


H = 0.3 


^l = 1.28(7) 


/ic — 2.6 


AtCT = 1.3 


H = 0.5 


^l = 1.52(7) 


Ate = 3.0 


Atcr = 1.5 


H = 0.8 


^l = 1.86(12) 


Ate — 3.6 


AtcT = 1.8 



TABLE I; Summary of the exponents found for the stress/contact ACFs at small /c for different roughness exponents H and the 
different methods analyzed in this work. Results for GW and Persson-^ are analytic predictions. GFMD results are numerical 
fits assuming ii^ = fj,c and the numbers in parentheses are uncertainties in the last significant digit. 



D. Comparison to Persson theory 

When comparing the GFMD results to Persson's theory, it is more convenient to compare the stress ACFs rather 
than the contact ACFs. Figures [6] and [7] show our numerical results for Ca-{q)) normalized by the stress ACF from the 
full contact approximation {f{q) = 1 in Eq. [7]). Results are shown for a wide range of loads that give relative contact 
area /c between 2% and 96%. Note that theory and simulation should not be compared for q > qs where C{q) ^ 
and the plotted ratio is not well defined. 
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FIG. 6: (Color online) The stress ACF, Ccr{q), normalized by the stress ACF for full contact as a function of q/qi for different 
relative contact areas Ac/Aq and roughness exponent H = 0.3. Results were averaged over the direction of q and over small 
bins in the magnitude to reduce numerical scatter. From Eq. [7] this ratio should be f{q), the relative contact area at a given 
magnification q/qi- Dashed lines show a numerical evaluation of this quantity for the same surfaces. 

Our numerical results for /c > 0.9 are nearly indistinguishable from the full contact expression. As fc decreases, 
the discrepancies from full contact increase. The magnitude of Ca{q) is depressed and the data drop with an apparent 
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power law indicating that /i > 2H. The deviation is clearly smaller for H = 0.8 than H = 0.3, but is present in both 
cases. 

In Persson's theory with corrections for partial contact (Eq. [7|) the ratio plotted in Figs. [6] and [7| should equal f{q), 
the relative contact area for surfaces resolved to q Dashed lines in the figures plot the numerically determined 
f{q) for the same surfaces. For smaller relative contact area, the theory still captures the trend correctly, i.e., the 
power laws with which the ACFs decay at large q are very similar to the numerical data. The theoretical prefactors are 
generally better for large H than for small H, which is consistent with the observation that the value for k predicted 
by Persson improves with increasing 
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FIG. 7: (Color online.) Same as previous figure, except that now the roughness exponent is H — 0.8 instead of H — 0.3. 
Note that all configurations were produced with the same random seed, which explains the correlation in the noise for diflterent 
Ac/Ao. 

For a few of the largest area fractions (not included here to make the figures clear) , there appears to be a crossover 
at a wave vector q* . For q < q* the results converge to the full contact prediction, while for q > q* the results follow 
the larger exponent observed at small area fractions. It is intriguing to note that the crossover behavior is observed 
in the simulations only near and above the percolation probability for a square lattice fc ~ 0.59^^, This suggests that 
when the contacting area percolates, the system behaves as if it were in full contact at large scales. The wavelength 
corresponding to q* might then correspond to the size of the largest non-contacting regions, leading to an increase in 
q* as the area fraction increases. Tests of these conjectures are beyond the scope of this work. 



E. Comparison to field-theoretical approach 

In referencei^, it was argued by one of us (MHM) that Persson's theory corresponds to the leading-order term 
of a rigorous field-theoretical expansion. The expansion is formally based on the assumption that a (free) energy 
functional exists describing the interaction between two contacting solids, which depends on the gap separating the 
two solids. 

For exponential repulsion, corrections to Persson's theory were worked out explicitly up to harmonic order. The 
main result relevant for the calculation of the correlation function is that equation (O, which is valid in Persson's 
theory, will be replaced with 

(a*(q)a(q)) ^ {Gi(q)}' ^q^ {h* {cDhid)) . (20) 

Here, G'i(q) = 1/(1 + CqE' /2(7o) is a correction factor that depends on the characteristic screening length ^ of the 
exponential repulsion, the magnitude q of the wavevector, the effective elastic constant E' and the macroscopic normal 
stress (To. As the normal force is never exactly zero, we have used f{q) = 1 for all values of q in the evaluation of 
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Eq. ([7]). Thus, like Persson theory, the field theory has no adjustable coefiicient. In the limit C ^ 0, corrections 
disappear. 

In figure [8] comparison is made between numerical results and the field theory. The numerical data was based on 
the same calculations as those presented in reference^-. It can be seen that corrections to Persson theory are very 
small for the smallest value of the screening length analyzed here, i.e., for ( = 0.01. However, for a value of C = 0.25 
the agreement between predicted and calculated stress ACF is essentially perfect for the relevant wave vectors. The 
degree of agreement is surprisingly good, given the relatively poor agreement in the stress histograms for that same 
value, see figure 1 in reference^^. At the largest value of (, that is, C = 1; there is almost perfect agreement, as to be 
expected for a rigorous expansion that is most accurate for large values of C- 




q/q. 



FIG. 8: (Color online) Calculation of stress ACFs for exponentially repulsive walls. Symbols indicate simulation data for 
different values of the screening length as measured in the units of the distance a between two discretization points. The 
dashed line reflects Persson's theory. Solid lines are drawn according to the field-theoretical approach presented in reference^^. 
Neither theory has adjustable coefficients. Here Ai = 1024o, As = 32a, H = 0.3, and L/Ao — 0.004. 



F. Analysis of connected contact patches 

A merit of GW and related theories is that they provide an intuitive explanation for why total load L and true 
area of contact are proportional to one another at small loads and thus that is constant. Although the relation 

2/3 

between area Uc and local load Ic for any individual contact is non-linear in these theories, Uc oc Ic , the number of 
contacts of each area rises linearly with load. As the total load increases, a contact that already exists will have an 
increased local load and grow in size. However, new contacts will form under increasing L so that there is a supply of 
new contact patches with small local loads. Under certain favorable conditions, the distribution of contact loads and 
sizes maintains the same shape and Uc = L/Ac remains constant. Previous work shows that while the distribution 
of contact areas is different than the GW prediction, it is independent of load at small loads®. Here we examine the 
relation between load and area within these patches. 

In figure[5]we present data for a number of systems with H = 0.3. In particular we analyze the load Ic that connected 
contact patches carry as a function of their area Oc- It can be seen that the data decomposes into two scaling regimes. 
At small Oc, the data follow the prediction of Hertzian contact mechanics that is used in GW, Oc oc Ic^^. At large 
flc, the load exhibits the linear scaling with area that is found for the entire macroscopic contact area. The crossover 
occurs when the contact area is comparable to the square of the small wavelength cutoff in the fractal scaling. At 
smaller scales single asperities are fully resolved, while at larger scales one sees the effects of roughness with many 
wavelengths. GW theory does not include the effect of these larger scales and it is interesting that the linear scaling of 
area and load enters so close to Ag. Hyun and Robbins have shown that there is also a crossover in the the probability 
of finding clusters of a given size at Uc '■-^ Xg. The probability is nearly flat for Uc < A^ and falls off as a power law at 
larger Uc ■ 
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FIG. 9: Load Ic that an individual connected contact patch carries as a function of its microscopic surface area Oc. Here 
H — 0.3, Ai — 64As = 4096a and the area fraction is indicated in the figure. 



IV. DISCUSSION AND CONCLUSIONS 



In this work, we computed the stress and contact ACFs that one obtains when pressing an elasticaUy deformable soHd 
against a rough, rigid, non-adhesive, and impenetrable substrate. Analytic arguments were presented for approximate 
bounds and exact sum rules on the stress and contact ACFs. These imply that decays more rapidly with wavevector 
q than Cc, but that the change in their ratio is only about a factor of two, no matter how large the system. As a 
result, the decrease of these ACFs with wavevector must be described by the same scaling exponent /u. 

Numerical results for /i were compared to the predictions of analytic theories for roughness exponents 0.3 < H < 0.8 
that span the typical range for experimental surfaces (Table|l|. Our GFMD results are consistent with ^ = 1+H within 
an error of 0.1. This is only half the value predicted by the bearing area and GW models, which neglect the elastic 
interactions between deforming asperities. Persson's theory includes these correlations through an approximation 
for the stress ACF that becomes exact in the limit of full contact. Recent extensions of his model^i^^ that include 
corrections for partial contact lead to a scaling exponent that is consistent with our numerical data22,. The prefactor 
predicted by this model appears to be slightly too high at small H , but to approach the numerical results as 77 — > 1. 
This is also the limit where the assumption of full contact is most accurate. In this context it is interesting to note 
that the value of k (Eq. [T]) also seemed to approach Persson's prediction as — > 1 in earlier work*. 

In this work we also tested whether a recently suggested field-theoretical approach to contact mechanics allows one 
to improve predictions for the stress ACFs^^. The new approach can be interpreted as an expansion, in which the 
perturbation parameter is a screening length which describes the exponential repulsion between two surfaces. A zero 
screening length corresponds to hard wall interactions. In this limit the leading order term of the expansion reduces 
to Persson's theory. Our numerical results suggest that including the next non-vanishing term vastly improves the 
agreement between calculated and predicted stress ACFs. 

Lastly, analysis of the load that is carried by connected contact patches revealed a crossover at a critical patch size. 
Smaller contacts exhibit a Hertzian relation between area and load. Larger contacts exhibit a linear relation between 
area and load. This linearity at the contact scale may be part of the reason that a linear relation between area and 
load is observed for the entire surface. 
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Here, we wish to remind the reader that the term "numerically-exact method" stands for a numerical method, in which a 
model, such as an elastic manifold pressed against a corrugated wall with hard-wall repulsion, is solved in such a way that 
all systematic errors can be controlled, e.g., by the number of mesh points. In contrast, the GW and Persson approaches 
make uncontrolled approximations. 

Note that a is calculated from the force from the rigid substrate normalized by the area, a^, per "atom" in the GFMD. Our 
calculation actually uses the z component of this force which differs from the normal component by a factor of the cosine of 
the angle between the surface normal and the z-axis. Since the rms slope is only 0.031, the difference is negligible. 
We found fic was 3.1 for Ai/As ~ 64 and 2.8 for Ai/As = 2048 for the bearing area model, while the size effect was merely 
0.1 for the GFMD data. 
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